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The elastic moduli of vortex crystals in anisotropic superconductors are frequently involved in the 
investigation of their phase diagram and transport properties. We provide a detailed analysis of the 
harmonic eigenvalues (normal modes) of the vortex lattice for general values of the magnetic field 
strength, going beyond the elastic continuum regime. The detailed behavior of these wavevector- 
dependent eigenvalues within the Brillouin zone (BZ), is compared with several frequently used 
approximations that we also recalculate. Throughout the BZ, transverse modes are less costly than 
their longitudinal counterparts, and there is an angular dependence which becomes more marked 
close to the zone boundary. Based on these results, we propose an analytic correction to the nonlocal 
continuum formulas which fits quite well the numerical behavior of the eigenvalues in the London 
regime. We use this approximate expression to calculate thermal fluctuations and the full melting 
line (according to Lindeman's criterion) for various values of the anisotropy parameter. 



I. INTRODUCTION 

Among the fascinating aspects of the new high- 
temperature superconducting materials, there are many 
which aije due to the rich and complex behavior of vor- 
tex linesEl. These quanta of magnetic flux penetrate the 
superconductor above a certain threshold value of the ex- 
ternal magnetic field H, the so-called lower critical field 
i?ci (roughly lO^^T , and their concentration increases 
with H up to the upper critical field (anncoximately 
lO^T), above which the material is normaflcl (see Fig- 
ure 0). In a classic work, Abrikosov showed that the min- 
imum energy arrangement of flux lines in a conventional 
superconductor is a triangular lattice, with a latticejSpac- 
ing a which varies with the magnetic field strengtha. 

Recently, with the discovery of high temperature 
cuprate superconductors, there has been a resurgence of 
interest in the properties of vortex matter. One novel 
aspect of these materials is the symmetry of the super^. 
conducting order parameter. A series of experimentsEl 
seem to provide evidence consistent with d-wave symme- 
try, in contrast to conventional superconductors which 
have s-wave symmetry. In addition to modifying the in- 
ternal structure of the vortex lines, d-wave symmetry has 
implications far the global mean-field arrangement of the 
vortex latticed. In particular, an oblique square lattice 
appears to replace the standard triangular arrangement, 
as the most stable configuration in part of the B — T 
phase diagramQ. (The new phase appears for high values 
of the magnetic field, whereas the triangular arrangement 
is confined to low fields.) 

Furthermore, the relative flexibility of vortex lines in 
these materials makes them susceptible to distortions by 
thermal fluctuations, and other sources of disorder (oxy- 
gen impurities, grain boundaries, etc.). The regular lat- 



tices obtained in mean-field theory are thus distortjei 
giving rise to a rich variety of vortex-matter phasesu'l 
It is thus necessary to understand the elastic response 
of vortex lattices to distortions, a subject that has been 
intensely studied in the context of conventional super- 
conductors. The stability of the triangular lattice against 
small distortions is guaranteed as long as the characteris- 
tic length of the field variations, the so-called penetrati^ 
depth A, remains smaller than the size of the systemll3. 
(An infinite A renders the lattice of vortex lines unstable 
against fluctuations or, more precisely, against shear de- 
formations.) At low temperatures, fluctuations are well 
described by small corrections to the mean-field results. 
Although the phenomenology is similar for d-wave su- 
perconductors, less is known about the effects of thermal 
fluctuations and disorder in this case. We shall thus fo- 
cus on a triangular lattice of vortex lines. The extension 
of the results to oblique configurations is possible, and 
should prove interesting. 

The elastic properties of the triangular vortex lattice 
at long distances are characterized by its compressional, 
shear, and tilt moduli. These moduli are frequently in- 
volved in the theoretical and experimental determination 
of the properties of the material, as for example, the in- 
tricate vortex phase diagrams. As such, the derivation 
of the elastic energy and the elastic jXaoduli has been un- 
dertaken in several publicationalJtEl. While it is well 
known that the elastic moduli of the vortex lattice de- 
pend strongly on the magnetic field strength, the values 
currently presented in the literature are only strictly use- 
ful in certain limiting situations. Most frequently, results 
are obtained in the so-called continuum limit, in which 
one disregards the discrete nature of the underlying vor- 
tex lattice. Obviously, this description is not suitable 
for the whole range of possible flux-line densities in the 
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mixed state, which extends from _ffci up to the upper 
critical field iJca , at which the magnetic field penetrates 
uniformly into the material. The elastic properties of 
the vortex lattice, and hence its stability against ther- 
mal fiuctuations, depend sensitively on the value of the 
magnetic field in this range. Furthermore, the important 
fluctuations sometimes occur at short wavelengths, where 
a simple elastic description may not be appropriate. It is 
thus worthwhile to obtain the general dependence of the 
energy cost of harmonic distortions of the vortex lattice. 
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FIG. 1. Schematic mean-field phase diagram of a type-II 
superconductor. The numbers are indicative orders of mag- 
nitude for highly anisotropic materials like YBa2Cu307 and 
BiaSraCaCuaOg. 

As temperature is increased, thermal fluctuations 
cause the vortex lattice to melt into a vortex liquid. 
Several experiments employing quite different techuiqi. 
have provided firm evidence for such a transitioncill 
In the new cuprate superconductors, the melting transi- 
tion can occur at temperatures well below the mean field 
point, so that the Abrikosov lattice is melted over a sub- 
stantial portion of the phase diagram. Furthermore, the 
vortex lattice can melt not only by increasing temper- 
ature, but also by decreasing the magnetic field to the 
vicinity of Hc^{T). In this region of the phase diagram, 
the concentration of vortices is very dilute; the separation 
a between neighboring flux lines is larger than the pene- 
tration length A, and the vortex- vortex interaction decays 
exponentially. As a consequence, the elastic moduli be- 
come exponentially small, and correspondingly, thermal 
fluctuations are greatly enhanced. This behavior gives 
rise to an interesting reentrant behavior of the melting 
line at low flelds. 

In this paper, we take up the computation of the elas- 
tic properties of a vortex lattice, in a systematic manner 
which is not restricted to the most commonly considered 
continuum limit. Our analysis allows the illustration of 
some of the subtleties involved in the calculation, and to 
successfully resolve some of the unclear points that we en- 
countered in our review of the literature, and which may 



also have puzzled other investigators of the subject. The 
paper is organized as follows: In Sec. |l|, we introduce the 
Ginzburg-Landau Hamiltonian for an anisotropic super- 
conductor. This model is commonly used in the literature 
as the starting point for studying the effects of fluctua- 
tions and disorder. Minimization of this Hamiltonian 
gives the optimal lattice conflguration, around which we 
then study the cost of distortions in the harmonic ap- 
proximation. This section also introduces the main pa- 
rameters and notation used throughout the paper. In 
Sec. [II, we provide an expression for the harmonic ker- 



nel in terms of a sum over Bravais lattice vectors or, 
equivalently, over reciprocal lattice vectors. The nor- 
mal mode eigenvalues are explicitly written in terms of 
this kernel. Sec. |^ is composed of subsections in which 
we introduce certain limiting situations, corresponding 
to a single line, local, local-continuum, and nonlocal- 
continuum limits. These limiting cases are often quoted 
in the literature, and are widely used in studies of vor- 
tex matter in high temperature superconductors. Within 
our framework, we recalculate the analytic values of the 
elastic eigenvalues in these limits before discussing our 
more general results in Section ^ where the harmonic 
eigenvalues are numerically evaluated. Our results are 
graphically presented in several plots, and compared to 
the limits introduced in Sec. |v| , to easily visualize the 
accuracy of the approximations involved. As a practical 
illustration of the potential applications of our results, 
we use the harmonic eigenvalues to calculate the ther- 
mal distortions of the vortex lattice in Section 0. The 
leading contribution to flux-line fluctuations in real space 
comes from the transverse modes. In conjunction with 
the Lindcman's criterion, we can then flnd the full form of 
the melting line as a function of the magnetic field. This 
is one of many potential applications that our general 
analysis makes possible, without the need to extrapolate 
the elastic behavior of the lattice to regimes beyond their 
range of validity. Finally, in the last section we summa- 
rize our main conclusions. 



II. GENERAL FORMULATION 

Our starting point is the continuum Ginzburg-Landau 
free energy for an anisotropic superconductor in the Lon- 
don limita. In this limit, the penetration length A (~ 10'^ 
A) is much larger than the coherence length of the super- 
conductor ^ (~ 10 A), and fiuctuations in the magnitude 
of the order parameter '^o are neglected. The phase de- 
gree of freedom is then the only variable that needs to 
be considered. For the anisotropic superconductors un- 
der consideration, this approximation breaks down in a 
narrow band close to , where the separation between 
vortices becomes comparable to ^. The Ginzburg-Landau 
Hamiltonian in this limit is given by 
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Here 0(rj^, z) is the phase of the complex order param- 
eter field ^'(rj^,z) = \E'oe**'''^'^\ A(rj_,z) is the mag- 
netic vector potential related to the magnetic induction 
through b(rj^,z) = V x A(rj_,z), a = [Ti'^ oY /"^Tn, 
4>o = hc/2e is the flux quantum, and = ra/M is the 
usual anisotropy parameter, with m and M being the ef- 
fective masses in the Cu-0 planes of the material, and 
along the perpendicular c axis chosen to coincide with 
the z direction, respectively. The external magnetic field 
H is also oriented parallel to the z axis H = Hiz ■ Note 
that Vj_, r±, and Aj_, denote the planar components of 
V, r, and A, respectively. 

Consider a configuration with N vortex lines, and de- 
compose the corresponding phase field into two parts: 
S^=i^n['''-L — Rn,±{z), z], which represents the singu- 
lar phase due to the N vortices passing through points 
Rn{z) = Rn,±{z) + zz {n = 1, N) on independent 
planes at each z; and 9^, a regular phase field account- 
ing for the relaxation due to the couplings between the 
planes. By construction, each 9!^ is the solution of a 
two-dimensional problem with the circulation constraint 

d9^ = 27r on any closed circuit C around the n-th vor- 
tex. We have chosen the coordinate z to parametrize the 
trajectories of the different flux lines. Nonetheless, one 
has to bear in mind that the results should be invariant 
under an arbitrary reparametrization. 

Variations of Eq. (|l|) with respect to the phase 9^ , and 
the vector potential A, provide the differential equations 
for these quantities, whose solutions minimize the energy 
in After considering the Coulomb gauge V • A = 0, 
these equations read 
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where k = q + k^z, Ac — A/e is the penetration length 
along the c-axis of the superconductor, and we have in- 
troduced the function 



2m 

r,1 



(,x.)./..e*=e--.^<- 



dR.a.x{z) 



dz 



(7) 



which represent the Fourier transform of the vortex func- 
tion dz9!^. 

Because of its singular nature, the Fourier transform 
of V^^^j, actually has a component 
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which is perpendicular to q. We thus define the transver- 
sal and longitudinal components of the in-plane projec- 
tion of the gauge field, as = (/ — qq) ■ Aj_, and 
Aj^ = qq ■ A± respectively. The former is obtained from 
Eq. (I) as 
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while the latter is obtained from the Coulomb gauge con- 
dition as 
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For a given distortion of the flux-lattice, the above so- 
lutions provide the form of the phase and gauge fields 
that minimize Eq. (0). The next step is to evaluate the 
energy cost of such distortions by substituting these so- 
lutions into the Ginzburg-Landau free energy, resulting 
in 
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where we have introduced the squared penetration 
depth in a plane perpendicular to the z axis, A2 = 
(/)2/(327r3a) = mc2/(167re2^'2), and A_l stands for the 
in-plane component of the Laplacian operator. The cou- 
pled equations (||) and (H) are easily solved in Fourier 
space to give 
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Let us now assume that the position vector Rn.i_ {z) is the 
sum of a perfect lattice vector, plus a small displacement 
field due to fluctuations, i.e. Rn,jL{z) = ii° -|-m„(z). Up 
to second order in the displacements u„(z), the energy 
cost can then be expressed as 



T-L^l-L° + AH. 



(12) 



The first term H°, gives the free energy of an array of N 
straight flux lines oriented parallel to the external field, 
and located at positions 
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Here L is the sample thickness, Hc^ — 0o/(47rA2) ln(K) 
is the lower critical field, k = A/^ is the Ginzburg- 
Landau parameter, H is the magnetic field strength, 
to — '/'o/(47''A)2 is an interaction energy per unit length, 
Ko{x) is the modified Bessel function of zeroth order, 
and dnm = Rn ^ R-m S'''*^ ^hc relative position vectors 
of any pair of flux lines. On a perfect triangular lat- 
tice, these vectors are of the form dnm = a{nei + me2), 
(n, m = 0, ±1, ±2, . . .), with a the lattice spacing, ei ori- 
ented, for instance, along the x axis, ei = i, and 62 = 
cos(7r/3)i;-t-sin(7r/3)y (see Fig. 0). The modulus of one of 



these vectors is then given by dnm — a{n 



nm) 



1/2. 

The first term on the right hand side of Eq. (13) rep- 
resents the energy cost of a single vortex line in a type 11 
superconductor, times the number of lines N. The last 
term is due to the interactions among fiux lines, which 
naturally depend on the interline separations. The pene- 
tration length A, sets the extent of the interaction poten- 
tial Ko{x), which diverges logarithmically at short dis- 
tances, and decays exponentially for a; 1. 

The quantity ATi. represents the harmonic contribu- 
tion of fiuctuations to the free energy. After writing the 
displacement fields Un{z) in terms of Fourier modes 
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and taking advantage of the translational symmetry of 
the lattice, the energy AH can be written as 
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where the index BZ indicates that the integration is per- 
formed over the first Brillouin zone in reciprocal space. 

All the relevant information is therefore contained in 
the harmonic kernel Mai3{Q,kz). Formally, we have 
to find the eigenvalues and eigenvectors of this matrix, 
and we can then calculate the extent of fluctuations in 
real space, the fluctuation corrections to the free en- 
ergy, and other relevant quantities. In Fourier space, the 
eigenvectors of M^piQ^k^) are N longitudinal modes, 
u^{Q,kz) = Q • u{Q,kz), and N transversal modes, 
u^{Q,kz) = IQ X u{Q,kz)\, with corresponding eigen- 
values h-hiQ^kz) and h-T{Q,kz)- In practice, it turns 
out that the analytic expressions for these eigenvalues is 
rather complex. That is why in the literature only certain 
limiting regimes are usually treated. We shall describe 
these limits later on in Sec. IV, and then go on to discuss 
their accuracy in comparison to our more general results. 



III. LONGITUDINAL & TRANSVERSAL MODES 

The longitudinal elastic eigenvalue A^iQ^kz) is typ- 
ically expressed in the literature in terms of the so- 
called compressional, cii{Q, kz), and tilt, C4^i{Q, kz), elas- 
tic moduli as cii{Q,kz)Q^ + CA4^(Q,kz)k^. (In general, 
however, we shall see that this decomposition may be 
inaccurate.) In terms of the matrix Mq,^(Q, /c^), the lon- 
gitudinal eigenvalue is given by 



AL{Q,kz) = QaMo,f3{Q,kz)Q0, 



(16) 



where, as usual, a repeated index is summed over. On the 
other hand, the transversal eigenvalue At{Q, kz), com- 
monly written in the literature in terms of the shear, 
ceeiQ, kz), and tilt moduh as CQe{Q, kz)Q'^ + cn{Q, kz)kl, 



At{Q, kz) = M^,,[Q, kz) - Qc.Mo.ii{Q, kz)Qp. (17) 

There are two alternative ways of expressing the in- 
teraction kernel Map(Q, kz): in terms of a sum over the 
Bravais lattice vectors, or as a sum over the reciprocal 
lattice vectors. The former yields 
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where we have introduced the areal density of flux lines 
n = N/A. For compactness of notation, let us also in- 
troduce the dimensionless variables x — -^/A^fc^ + 1, and 
Xc = e X. In Eq. (ITsI), the quantity neo £{kz)/X^ repre- 
sents the line tension of each vortex, which, in general, 
contains contributions from both the Josephson coupling 
between the different Cu-0 layers in the material, and 
the magnetic interlayer interactions, as 
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(19) 



The interaction kernel Rafi{d, kz) is given by 
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Note that we use the symbol d — \d\ to indicate the 
distance between a pair of lattice points. The sum in 
Eq. (|l^) runs over all such separations from a specific 
point on the lattice, as illustrated in Fig. ||. 
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FIG. 2. Lattice vectors with a common origin at the lattice 
point P, and with the same modulus d. 

The matrix R{d, kz) depends on the different relative 
position vectors d = dnm of the perfect triangular lattice. 
When summing over all lattice vectors with a common 
origin at a point P, we can take advantage of the lattice 
symmetries. For example, the sum over products of an 
odd number of vector components vanishes, as in 



^ da = 0, 
d#0 



and. 



d#0 



dadpd^ — 0. 



(21) 



On the other hand, the sums involving products of an 
even number of components are nonvanishing and con- 
strained by the symmetries; for example. 



^dadp = 5af3^ -77, 



(22) 



d#0 



dadpd^ds — [SafsSyS + Sa-fSpS + SasSpy) ^ — . 

(Higher order terms have a more complex structure.) The 
sums J2d^o right hand side of Eqs. (22), are equiv- 

alent to J2d^o with a degeneracy factor g{d), which 
counts the number of vectors whose modulus is d. 

Equivalently, we can express the kernel as a sum 
over reciprocal lattice vectors G. The same symme- 
tries (pl|)-(|2^) are, of course, valid for these vectors, 
which, according to our choice of dnm, are given by 
G = 47^/(^/3a)(p5l + gga), (p, q = 0, ±1, ±2, . . .), where 
47r/V3 a is the reciprocal lattice spacing, and the unit 
vectors are gi = sin(7r/3)i — cos(7r/3)?/, and 32 = y- In 
terms of these vectors, the interaction kernel reads 

2 J 2 

M^piQ, kz) = ^[Gf^(G - Q, kz) ~ G^p{G, 0)] 

^ G 

j:^[Gip{q-QM-Gc.p{q,kz)]]- (23) 

The symbol ^ indicates symmetrization with respect to 
Q, i.e., Gip{G - Q, kz) = [G„^(G - Q, kz) + G„^(G + 
Q, kz)]/'2.. (The second term on the right-hand side of the 
expression is included to properly account for the partic- 
ular case n = m and 2: = z' in Eq. (pi]).) In addition, 
from Eq. (11) we obtain 



Our goal is to provide accurate values of the harmonic 
eigenvalues as a function of Q and kz within the BZ, 
and for different concentrations of flux lines. The di- 
mensionless parameter a/ A will be used as the indica- 
tor of the areal density of flux lines. A small value of 
a/ A, corresponds to a highly dense regime with strong 
and non-local interactions among the flux lines. On the 
other hand, for a/A ^ 1 the concentration of flux lines 
is very low, and interactions among them very weak. In 
this dilute limit, the elastic behavior of the lattice re- 
flects to properties of a single flux line. We thus evaluate 
numerically the harmonic eigenvalues as a function a/A 
throughout the BZ. The equivalence of expressions (|lq ) 
and (|2^) renders the use of one or the other a matter of 
convenience. For instance, for very high areal densities of 
vortices n, it is common to disregard the discreteness of 
the underlying arrangement and approximate the sums 
over lattice positions by integrals. This so-called-contin- 
uum limit is very often used in the literatureE3~E3. The 
elastic moduli which follow from this approximation can 
be read directly from Eq. (p4|), when only the recipro- 
cal lattice vector G = is taken into account. We shall 
comment on the accuracy of this limit in the following 
sections. 



IV. LIMITING REGIMES 

In this section, we first present some special cases pre- 
viously discussed in the literature. By comparing these 
limits to our numerical results, one can then see their 
range of validity and the accuracy of the approximations 
involved in their formulation. We shall also emphasize 
the roles played by the interline distance, the penetration 
length A, and the symmetries of the triangular lattice. 



A. Single line 

At very low areal densities, i.e. a/ A ^ 1, the kernel 
Maf}{Q,kz) should be just N times the result obtained 
for a single flux line. All the modifled Bessel functions 
appearing in Eq. (|20| ) decay exponentially fast for large 
values of their argument (which is proportional to a/ A), 
and only weakly contribute to the flnal values of the har- 
monic eigenvalues. In other words, in Eq. (18) only the 
line tension part £{kz) is important for large values of 
a/ A. In Sec. we show that the single-line approxima- 
tion indeed provides the best description of the elastic 
properties of the system at very low concentrations of 
flux lines. 



B. Elastic moduli 

The energy cost of long wavelength distortions is de- 
scribed by an elastic theory, whose form is constrained 
by symmetries. The triangular lattice is isotropic, and 
governed by the compressional, cn, shear, cge, and tilt, 
C44, moduli. In terms of the Fourier modes, the elastic 
energy is 
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+ cee{Q^Sai3 ~ QaQis)] Ua{Q, k;,)up{-Q, -fc^), (25) 
and the harmonic eigenvalues have the simple forms 

J^riQ. k.) = [ceeQ^ + Ci^kl)/2. (26) 



Formally, Eq. (gj) together with (g^ or (|2^) provide 
the harmonic eigenvalues throughout the Brillouin zone. 
Naturally, the elastic limit is regained by expanding these 
results up to second order in Q and kz. Neglecting the 
higher order terms is referred to as the local limit, as it 
is usually obtained by including only short-range inter- 
actions. Sometimes, non-local elastic moduli are intro- 
duced which depend on Q. This is not always useful, as 
it constrains the form of the harmonic eigenvalues as in 
Eqs.(p6|), whereas symmetry allows higher order powers 
of Q to appear in other forms. 

After expanding the cosine, and using the symmetry 
properties of a triangular lattice in Eqs. (|2l|)-(|2^), as 
well as cettain relationship among the modified Bessel 
functionstHI, we arrive at 
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where 8^{kz) is the single line tension in the local ap- 
proximation, £\kz) = (eAfc^)^ ln(K/e) -I- {Xk^Y /2 (valid 
for K > 1 > e). 

Comparing this expression with Eq.(p5|) allows us to 
identify the compressional modulus 
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the shear modulus 
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and the tilt modulus 
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Note that the mean-field lattice spacing a, is related to 
the strength of the applied field H through the relation- 
ship 



H = Hr 



47rA2 



-y 



Ko 



(34) 



which is obtained from dH°/da\a = 0. Then, by com- 
parison, the local tilt modulus for an isotropic material 
(e — 1) can be written as C44 — n(j)oH/'i'K^ in agreement 
with the result expected for the local^ilt modulus of a 
rotationally invariant superconductorlij. Strictly speak- 
ing, in Eq. ( |3^ ) there is an extra factor of 1/2, coming 
from thi; local line tension. However, as pointed au± by 
Fisherll3, one could chose a short distance cutoff in- 
side ln(A/f_) to reproduce exactly the local isotropic limit. 

In Figs, ga) and ||b), we depict the functions F^ and Fg, 
for different values of the dimensionless quantity d/\, to 
illustrate some subtleties associated with the estimation 
of the shear modulus, often forgotten or misunderstood 
in the literature. As shown in Fig. ^3), for small values 
of d/\ the function Fg is negative, whereas if d/\ > 2 it 
becomes positive and then decays exponentially to zero. 
The implications of this functional form for the shear 
modulus are as follows: As indicated before, the range 
of interactions among vortex lines is determined by the 
penetration length A. For a dilute vortex lattice whose 
lattice spacing a is comparable or even greater than A, 
almost all the terms in the sum in Eq. ( |30|) are positive, 
and as they rapidly decay to zero, only the first few lat- 
tice vectors for neighboring sites are needed to calculate 
the shear modulus. On the other hand, in a dense system 
there are many neighboring vortices within the interac- 
tion range A, contributing a negative amount to the sum 
in Eq. (pQ). It is then clearly not sufficient to consider 
only interactions among a few neighboring lines, and to 
account for the stability of a dense lattice against shear 
deformations, we need to sum over many lattice vectors 
(finally giving rise to a positive shear modulus). 
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FIG. 3. a) Behavior of Fc in Eq. (^9|) as a function of the dimensionless parameter d/X. b) Behavior of the function Fs in 
Eq. as a function of the dimensionless parameter d/X. The long-dashed hues represent the A — > oo hmit of these functions. 



In Figs. |a) and|b), we also indicate with a dashed 
Une the hmit A — > oo (for a finite value of d) . In this case 
there is a logarithmic interaction anaang the vortices, as 
in the two dimensional Coulomb gag23. Note that in this 
limit, all factors in the sum for the shear modulus are 
negative, and the lattice is unstable to shear deformation. 
(This is a manifestation of the more general instability 
of a system of point charges in the absence of external 
potentials 

Figure Q depicts the compressional, shear, and tilt 
moduli as a function of a/ A, obtained by numerical evalu- 
ation of the sums in Eqs.(p8|), (^o|), and (^^, respectively. 
These moduli have been normalized by their respective 
values in the local continuum limit, which for the long 
wavelength compressional modulus is defined as 



and, for the tilt modulus we obtain 



C44 



A 



-Ft 



pc 



4^A2 



A 



pc 



(37) 



Note that for the latter, the term d 
eluded in Eq. (||). 

Similar results for the continuum 1 
previously reported in the literature 



is already in- 



it have been 
In Ref. |l|. 



Cll 



neo 



-Fc{d ^ 0) 
167rA2 



-1 



B<t>o 
47r (87rA)2 ' 



(35) 



where we have replaced the sums by integrals, after first 
adding and susbtracting the d ^ element explicitly, 
and introducing the appropriate unit area Ape = V?ia? /2, 
i.e. the area of the primitive cell. We have also de- 
fined the average magnetic induction B = ncpo, with 
n = N/A = l/Apc- Similarly, the long wavelength shear 
modulus is 



nCo I / , s / d d ( d 



A^ 



Kogan and Campbell calculated the shear modulus of a 
flux lattice in different situations. Essentially, they pro- 
vided an expression for the shear modulus cee in terms of 
a sum over components of the reciprocal lattice vectors. 
Because of the symmetry relations of Eq. ( p^ ) quoted in 
the previous section, all the terms included in their sums 
cancel each other out, and should yield a shear modulus 
equal to zero. Kogan and Campbell numerically evalu- 
ated this sum over r ecip rocal lattice vectors, and reported 
a value close to Eq. (|3^) . They argued that there is a con- 
tribution-due to the specific form of the cut-off at short 
distancestj. However, in view of the above derivation of 
Eq. (p6[) , it is clear that cutoffs are not relevant to the 
local continuum value of the shear modulus. 

As soon as the lattice spacing is greater than the pen- 
etration length, a > A, all the compressional, shear, 
and tilt moduli deviate strongly from the continuum re- 
sults. The compressional and shear moduli then decay 
exponentially to zero, while the tilt modulus attains a 
constant value of neo(e^ ln(K/e) + 1/2). (This is why 
when normalized by its continuum value of Eq.(^7|) the 
ratio grows as (a/A)^.) As we shall demonstrate in 



!^(i + o) = -^ 

4 ' (87rA)2' 



(36) 



Sec. (VI), under these circumstances the lattice again be- 
comes rather soft and liable to melting, as its fluctuations 
are greatly enhanced for large interline separations. 
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FIG. 4. a) Normalized compressional and shear moduli as a function of the dimensionless parameter a/A, obtained by numerical 
evaluation of the sums in Eqs. (^) and (^o|). b) Normalized tilt modulus as a function of the dimensionless parameter a/A, 
from numerical evaluation of the sum in Eq.d32). 



C. (Non-local) Continuum limit 

The continuum limit refers to a situation in which the 
concentration of flux lines is high. In this limit, one disre- 
gards the underlying triangular lattice arrangement-jafld 
treats the array of vortices as a continuous mediurnEZrEj. 
In the continuum limit, the sum of Eq. ( p^ ) over Bra- 
vais lattice vectors is approximated by an integral, or 
equivalently, in the sum over reciprocal lattice vectors 
in Eq. (|2^) only the vector G = is considered. As a 
result, the elastic kernel reduces to 



Ma/3(Q,fc.) = 



(A2Q2 + A2fc2 + 1 

1 /■ d^q 



[X^Q^ + kl) + l] Q2 

e QciQf) 
>al3 



(38) 



The last term in this equation is a rather intricate inte- 
gral, which, on the other hand, can be easily evaluated in 
the local approximation, i.e. when keeping terms up to 
and only. With this approximation, we can rewrite 
Eq. (H) as 



Ma(3{Q,k^ 



(Q2 + fc2) Q^Q^ 



kl 



[A2(Q2+fc2) 
QaQp 



1] 



{XIQ-^ + A2fc2 + 1) 

1 Q^Sal3 QcxQp 

n 167rA2 ~ 87rA2 



'a/3 



(39) 



By comparing with Eq. (p3) , we can now directly read off 
the nonlocal continuum (again, up to the approximation 
carried out to evaluate the integral) values of the elastic 
moduli as 



nl 
C44 



1 



47r (A2Q2 + A2fc2 + i)' 



(40) 



nl 
-11 



52 



1] 



47r [A2(Q2 + kl) + 1](A2Q2 + A2fc2 + 1) 



C66 



(87rA)2 



?^A)2' 
(41) 

(42) 



Note that, up to second order in and k^^ Eqs. (Efl) 



through ( |42| ) reproduce the local hmits calculated in the 
previous section from the real space expression of the in- 
teraction kernel. 



V. NUMERICAL RESULTS 

In this section, we present and discuss our results from 
numerical evaluations of Eqs. ([l6|)-(|l^) using, in particu- 
lar, the real space expression of the interaction matrix in 
Eq. (p^). In most cases, we selected a Ginzburg-Landau 
parameter n equal to 100, and an anisotropy ratio of 
e = 0.1. However, later in this section, we present some 
results for different values of e, which may be more ap- 
propriate to describe materials with stronger anisotropy. 
For the sake of simplicity, we introduce the dimensionless 
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eigenvalues = AL/ineo/X"^) and = At/ (neo/ X^), 
and discuss their behavior throughout this section. 




FIG. 5. Irreducible Brillouin zone of a regular vortex lattice. 



A. Angular dependence 

Because of the point group symmetries of the triangu- 
lar lattice, it is sufficient to consider only the irreducible 
Brillouin zone (IBZ), indicated in black in Fig. ||. In 
particular, we have chosen five paths corresponding to 
values of the angle a (defined within the figure) equal 
to 0, 7r/24, 7r/12, 7r/8, and 7r/6. Figure ^ shows the nu- 
merically evaluated eigenvalues along these paths within 
the IBZ. The horizontal axis is the dimensionless quan- 
tity {aQ)^. The plots also correspond to the choice of 
Xkz = 1.0. The dotted fines represent the longitudinal 
eigenvalue, A^, while we use dashed lines for the transver- 
sal eigenvalue, A^. 



The different plots in Fig. ^) clearly show that there 
is an angular dependence, which becomes more marked 
as we approach the edge of the BZ. This fact was already 
pointed out by Brandt in his, early work on the elastic 
properties of vortex latticeslla. The anisotropy is more 
pronounced in the dense regime, and diminishes as we 
increase a/A. In the concentrated case with a/A = 0.2, 
the relative differences |A(a) — A(7r/12)|/A(7r/12) are up 
to 20% for the longitudinal eigenvalue, and 60% for the 
transversal one. These ratios decrease in the intermedi- 
ate region with a/A = 1.0, and are of the order of 10% 
and 15%, respectively, for a/A = 5.0, within the dilute 
regime. We also observe that the longitudinal eigenvalue 
is largest along the a = 0-direction, while the transversal 
eigenvalue attains its smallest value for a = 0. 

Note tfiat both eigenvalues, for all angles a, have the 
same value at Q = 0, determined by Afc^ and a/A. (Since 
this is not immediately apparent from Fig. |^, the inter- 
section point is marked with a black dot on the vertical 
axis.) Away from Q = 0, the transversal eigenvalue drops 
sharply, more so in the dense regime. Transversal modes 
are therefore less costly than longitudinal modes. Fig- 
ure ^ also indicates that at a finite value of Afc^ (equal 
to 1.0 in this case), the minimum cost is obtained for a 
finite value of the rescaled wavevector aQ, which depends 
on the density of flux lines. 
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FIG. 6. Dimensionless elastic eigenvalues Ay and A^ as a function of (aQ) , for Afcz — 1.0, and at a) a/A — 0.2, b) a/A = 1.0, 
and c) a/A = 5.0. The dotted lines correspond to the longitudinal eigenvalue, while dashed lines indicate the transversal 
eigenvalue. The different lines are obtained for different paths along the IBZ, at angles a = 0, 7r/24, 7r/12, 7r/8, and 7r/6. The 
two solid lines in a) and b) depict eigenvalues from the nonlocal continuum limit in Sec. 

est concentration of flux lines (a/A = 0.2), and even at 
the intermediate concentrations of a/A = 1.0. However, 
the approximation fails at lower densities, and is not even 
included in Fig. ^c) , as the predicted eigenvalues fall well 
outside the range of our plot. 



The solid lines in Figs. |6^) and |^) are the results 
of the nonlo cal co ntinuum approximation introduced 
in subsection IV C|, i.e., A^ ~ (c"{Q^ 
At 



c2iki)/2 and 



(ceeQ^ + C44A:^)/2, with the appropriate param- 
eters. From these figures, one can conclude that (at least 
for Xkz = 1.0) the continuum approximation reproduces 
extremely well the behavior of both eigenvalues along the 
a = 7r/12 direction. This is certainly the case at the high- 
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B. Variations with 

To explore the fc^-dependence of the resuhs, in Fig. 
we plot the longitudinal and transversal eigenvalues as a 
function of {aQ)^ for a fixed angle a = 7r/12, and for a 
number of values of a/ A and Xkz- From these figures, we 
see that, not surprisingly, both eigenvalues increase as a 
function of Afc^ , and consequently, modes with higher kz 
are in general costly. For the smallest value of a/ A = 0.2 
in Fig. ^), the Xk^ dependence is most pronounced in a 
small region close to the BZ center. This region becomes 
larger as we increase either Xkz or a/ A. We also note that 
both eigenvalues show a linear dependence on (aQ)^ for 
a large fraction of the IBZ. The slope of this line is al- 
most constant within each plot, i.e. for a given value of 
a/A. The transversal eigenvalue for Xkz = is a straight 
line throughout the whole range of values of {aQY , as 
is the longitudinal eigenvalue for high enough values of 
Xkz- These features are qualitatively well accounted for 
by the continuum limit results. According to Eqs. (40), 
( pT| ) , and ( p^ ) , the transversal eigenvalue 



AT(g,fc.)^ -[566Q' + cJl(Q,fc.)fc^] 



Q2 + 2Tm^e^ 



A2fc2 



(A2Q2 + A2fc2 + 1) 



(43) 



reduces to cgeQ^/S for Xkz = 0, i.e., grows linearly with 
{aQY' , with a slope neo/Sa^. In addition, for suficicntly 



large values of Xkz » AcQ, c^i{Q,kz)kz reaches a satu- 
ration value, rendering the term cqqQ^ /2 as the leading 
contribution depending on Q. Similarly, the longitudinal 
eigenvalue 

Ai(g,fc,) ~ ]^[cf,{Q,kz)Q^ + ct{Q,kz)kl] 

" °[A2(Q2 + A;2) + 1] 8 ^ ^ > 

reaches a saturation value for Xkz 3> AQ, so that the 
term —ntoQ^/^ controls essentially the Q dependence of 
the eigenvalue at large values of Xkz- Eqs. (p3|)-(^) fit 
very well the numerical resuls obtained for a/A = 0.2 
and Xkz = 0.0, 1.0, 10.0 (we are not showing these results 
in Fig. 1?^) for the sake of clearity of the figure). For 
the larger value of Xkz — 50.0, although the qualitative 
form of the functions is quite well approximated by these 
equations, there are major quantitative differences with 
respect to the numerical results mainly due to a smaller 
value of both eigenvalues at Q = 0. In Fig. ^p) with 
a/ A — 1.0, there are important differences in the value 
at Q = already for Xkz = 10.0, as well as small devi- 
ations from the value of the slope of the linear part of 
the curves. These deficiencies of the continuum approx- 
imation are even more evident for a/A = 5.0 in Fig. |^). 
In the latter case, we also note the emergence of single- 
line characteristics: The eigenvalues are clearly mostly a 
function of kz, with only small variations with (aQ)^. 




{aQf [aQf {aQf 

FIG. 7. Elastic eigenvalues as a function of [aQ)'^ for a — 7r/12, and a) a/A = 0.2, b) a/A = 1.0, and c) a/A — 5.0. The 
different lines correspond to values of \kz = 0.0, 1.0, 10.0, and 50.0 only in plot a). Dotted and dashed lines correspond to: 
longitudinal and transversal eigenvalues respectively. 



To better emphasize this discussion, consider the elas- 
tic eigenvalues along the Q = 0-axis, i.e. the point high- 
lighted by a black dot in the previous figures, as a func- 
tion of Xkz- The results, for the same values a/X — 0.2, 
1.0, and 5.0, are plotted in Figs. ^) through c), and 
compared to the limiting cases described in the previous 
section. In particular, the dot-dashed line in these figures 



is the outcome of a nonlocal continuum approximation, 
whereas the dotted and the dashed lines correspond to 
the local continuum and the local results defined in sub- 



section tVB. The long-dashed line is the graphical rep- 
resentation of Eq. (|l9[), i.e., of the single line behavior 
expected at very low densities of flux lines. 
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FIG. 8. Elastic eigenvalues, A_l(0, fcz) = AT{0,kz), as a function of Xkz for Q — 0, and a) a/A — 0.2, b) a/A = 1.0, and 
c) a/A = 5.0. The numerical results indicated by solid lines are compared to the limiting cases in Sec. The dot-dashed, 
dotted, dashed, and long-dashed lines correspond to the nonlocal continuum, local continuum, local, and single line limits, 
respectively. 



C. Comparison to analytic approximations 

We observe in Fig. | that AHQ = 0, k^) = A^(Q = 

0, fcz) increase monotonously with Xkz, with the single- 
hne behavior of £{kz) in Eq. (|9|) eventually taking over 
for values of Afc^ greater than a characteristic Afc^, pro- 
portional to (ea/A)~^. The contribution of the nonzero 
reciprocal lattice vectors in Eq. (^3|) becomes relevant for 
Xkz > Xkz- In the nonlocal continuum approximation, 

1. e., considering only G = 0, A^Q = 0,kz) = A^((3 = 
0,kz) reaches a saturation value of 2ttii?€o at large values 
of Xkz- Indeed, for high areal densities of a/A = 0.2, our 
numerical outcome is very close to the nonlocal contin- 
uum approximation only up to Afc^- On the other hand, 
at low areal densities, such as for a/X — 5.0, the nu- 
merical data are obviously much closer to the single-line 
limit, than to any of the other forms, over the whole 
range of values of Xkz- As shown in Figs. |[;) and 0c), 
the variations of A^ and At with {aQY in the latter 
case are very weak. The eigenvalues throughout the IBZ 
essentially coincide with their values on the Q = 0-axis, 
which is set by the form of £{kz)- In the interval between 
these two values of a/A, Al{Q = 0,kz) — At{Q — 0,kz) 
gradually cross over from saturation type behavior to the 
parabolic- log dependence expressed in Eq. (^9|). The re- 
sults of the local approximation naturally fit well the be- 
havior of both eigenvalues for small values of Q and kz , 
as expected. 

In discussing Fig. 0, we noted that while the non-local 
continuum expressions capture the qualitative form of the 
numerical results, there are also important quantitative 
differences. Using the insights gained from the numerics, 
we shall now present an analytic form that corrects some 
of these discrepancies. One difference from the contin- 
uum results arises from the limiting slope in Fig. |7[ In 
the continuum limit, the term eventually deter- 

mines the slope of the linear regime of the transversal 



eigenvalues. However, as shown in Fig. |^), the actual 
shear modulus cqq can be quite different from the lim- 
iting value of C66 used in the continuum approximation. 
The differences in slope thus merely reflect the differences 
between cee an d Cee as a function of a/ A, as already dis- 
cussed in Sec. |VB. This deficiency of the continuum 
eigenvalues is thus removed by using the exact value of 
C66 from Fig. |a). 

A second difference is an overall shift of the eigenval- 
ues from the continuum limit prediction, which becomes 
more pronounced at larger Xkz- This clearly originates 
in the differences appearing already in Figs. ^) and b) 
for Al(Q = 0,kz) = At{Q = Q,kz) at Q = and suf- 
ficiently large kz- Since the continuum approximation 
uses only the term with G = in Eq. (p3|), we introduce 
a correction by replacing the sum over all the remaining 
reciprocal lattice vectors with an integral. These correc- 
tions result in the following expression for the transversal 
eigenvalue: 



1 B"^ 
AriQ.kz) ~ -[ceaQ'^ + cf^kl] + —{5o,fi - Qo^Qp) 

Z OTT 



BZ 



[G^p{G,kz)-G^p{GM- 



(45) 



In this equation, the correct local shear modulus, which 
depends on the field strength or lattice spacing a through 
Eq. (^o|) , is used in place of its continuum limit. Also the 
non-linear modulus C44 is corrected by the additional in- 
tegral over the nonzero reciprocal lattice vectors, prop- 
erly normalized by the BZ area Abz — 87r^/(\/3a^). To 
exclude the point at G = 0, the radial component of G in 
the above integral goes from C, defined by Abz — t^C^, 
to 00 (or, if necessary, to the short distance cutoff ^^^). 
After evaluating the integral, we obtain 

AT{Q,kz)-^[ceeQ^ + cf^kl] 
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e^/c? In 



In 



A2fc2 + A2C2 + r 



A2C2 



1 



(46) 



We have ascertained that this expression provides an 
excellent fit to the numerically obtained results for 
AriQjkz) along the a = 7r/12 direction. In the next 
section we provide explicit comparisons of the ana- 
lytic expression and numerics for different values of 
the anisotropy. The longitudinal eigenvalue along the 
a = 7r/12 direction can also be reasonably well fit- 
ted by a similar analytic expression, which replaces the 
term [ceeQ^ + c^lkl]/2 in Eq. ( ^6|) with the expression 
in Eq. (p4[), with a further substitution of for 
neoQ^/8 = ce6Q^/2, as 



Ai(Q,/c,) ~2™2g^ 



1 



°[A2(Q2 + fc2)^ 

neo2,2, / KVe2 + A2fc2 + i 

ne^ / A2fc2 + A"C^ + l 
'4A2 A2C2 + 1 



(47) 



D. Anisotropy 

We conclude this section by discussing the dependence 
of our results on the anisotropy of the superconductor. 
The results presented so far correspond to a particu- 
lar value of the anisotropy parameter, namely e = 0.1, 
which falls within the range (from 1/10 — L/5) reported 
for YBaaCuaOr (YBCO) in the literature^. However, 
smaller values of e, in the interval 1/100—1/50, character- 
ize highly anisotropic materials such as Bi2Sr2CaCu2 08 
(BiSCCO). For such highly anisotropic materials, the dis- 
creteness of their layered structure becomes important. 



and one may well question the validity of a three di- 
mensional Landau-Ginzburg description. An alternative 
model is a set of weakly (Josephson) coupled supercon- 
ducting layers. The applicability of the three dimensional 
description is typically assessed by comparing the coher- 
ence length along the c axis = witluthe distance d 
between the Cu-0 layers in the materialtl. A coherence 
length larger than the layer spacing usuallv-iustifies 
a continuous description along the c directionEj. This 
is certainly the case for YBCO, but not necessarily for 
BiSCCO. Nevertheless, extrapolating the results of the 
continuous description may also provide insights into the 
elastic properties of highly anisotropic superconductors 
such as BiSCCO. To this end, we compare results ob- 
tained for three different values of the anisotropy param- 
eter, namely e — 0.02,0.1, and 0.5, for the usual areal 
densities of a/A = 0.2, 1.0, 5.0. 

In Figs. ^ and |l^, we again plot the elastic eigenval- 
ues along the Q = 0-axis, and the transversal eigenvalue 
along the a = 7r/12 direction inside the IBZ, respectively, 
for different choices of e. There is clearly a strong depen- 
dence on e, mainly originating from the e-dependence of 
the line tension £{kz) in Eq. ([l9|). Naturally, the most 
pronounced tendency is the decrease in eigenvalues with 
e, reflecting the general softening of the vortex lattice. 
It is worth pointing out however, that the slopes of the 
linear part of the curves in Fig. ^ are independent of 
e, reflecting merely the fact that the shear modulus cge 
is independent of the anisotropy parameter. The dashed 
lines in the three graphs of Fig. |9| correspond to the an- 
alytic expression of Eq. ( p6| ) (for the appropriate param- 
eters a/ A and e), and clearly provide an excellent fit to 
the numerical data along the Q — axis. Similarly in 
Fig. 10, the solid lines are from Eq. ( ^6| ) for the transveral 
eigenvalue, again showing fairly good agreement with the 
numerical data. Analogous results are obtained for the 
longitudinal eigenvalue (not shown). 
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FIG. 9. Elastic eigenvalues along the Q = 0-axis as a function of Afcz, with different values, e = 0.02, 0.1, and 0.5, of the 
anisotropy parameter. The dashed lines depict the analytic expression in Eq. (Uq) with the corresponding choice of parameters. 
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FIG. 10. Transversal eigenvalue along the a = 7r/12 direction in the IBZ as a function of (aQ)^, and for Afcz — 10.0. The 
dashed lines are numerical results for e = 0.02 (BiSCCO), 0.1 (YBCO), and 0.5, while the solid lines depict the corresponding 
results from the analytic expression in Eq. (wyt). 



VI. THE MELTING LINE 

The melting of a vortex lattice by thermal fluctuations 
has attracted considerable attention in the context of 
high temperature superconductors. This transition has 
been observed by means of several experimental tech- 
niques such as bulk magnetization, local induction, and 
latent heat measurementsEiTcEl. The line-like nature of 
the constituent elements provides intriguing challenges 
to theoretical analysis. Important features of the melt- 
ing transition are the negative slope of the melting curve 
T,n{B) at high fields, its reentrant behavior at low fields, 
and its marked dependence on anisotropy. Some of these 
features can be extracted from sinmle models of the votj 
tex lattice, as in the so-called XiH, BosS^^ and cag^ 
models. These models agree in the prediction of certain 
universal features, such as the scaling of the melting tem- 
perature with the anisotropy parameter, or the ma^aetic 
field, in the high field region of the phase diagramL£l. 

In the absence of a rigurous theory for three- 
dimensional melting, the position and shape of the vor- 
tex lattice melting line is usually estimated using the so- 
called Lindeman criterion. According to this criterion, 
the lattice melts when thermally induced fluctuations of 
a lattice point become comparable to the lattice spacing. 
This condition can be written as 



1/2 



(48) 



where cl ~ 0.1 — 0.2 is the (empirically chosen) Linde- 
man parameter. The extent of fluctuations is measured 
by the autocorrelation function, {v?Lxq)). This quantity 
was calculated Jay Nelson and SeungEj iy. the local limit, 
and by Brandtllj and Houghton et aZ.li3 for the nonlo- 
cal continuum limit, and in more general circumstances 
including variations of the amplitude of the order param- 
eter in the Ginzburg-Landau hamiltonian. 

In the preceeding sections, we calculated the harmonic 
energy cost of fluctuations of the vortex lattice. In a clas- 



sical equilibrium state, each independent harmonic mode 
acquires a thermal energy of ksT /2. By adding the cor- 
responding squared amplitudes of the normal modes, the 
autocorrelation function is obtained as 
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kaT f dk 
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(49) 



The Lindeman criterion now provides a rough estimate 
of the melting temperature T^, through the relationship 



kflTrr 
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d^Q 
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+ A^i) 



2 2 



(50) 



The leading contribution to the fluctuations comes 
from the transversal modes, which as discussed in Sec. IV 
are always smaller than the longitudinal modes. We 
also ignore the angular dependence observed in Figs. ||, 
which is rather weak, and should not give rise to qualita- 
tively different results. Furthermore, we assume that the 
transversal eigenvalue is isotropic, and given by Eq. (|4^). 
This expression can be integrated analytically over the 
two-dimensional BZ, resulting in a rather long expres- 
sion, that we have finally integrated numerically over Xkz- 
The numerical integral over Afc^ is, in all cases, between a 
long wavelength cutoff equal to i£r^ and the e-dependent 
short scale cutoff of A/^c = K/a3. 

Provided that we are sufficiently far away from the 
critical temperature Tc, one can further assume that the 
penetration length is independent of temperature, i.e. 
A(T) ~ A(0). As we get closer to Tc, however, one must 
consider the temperature dependence of A in order to 
obtain sensible results. For the sake of consistency of 
our analysis, within the Ginzburg-Landau critical regime, 
we then use the mean-field temperature dependences of 
both A and ^, i.e. A(T) = A(0)(1 - T/Tc)-^/'^ and 
C(r) = ^(0)(1 - TjTc)-^!'^. We expect this approxima- 
tion to fail in close vicinity of due to critical fluctua- 
tions. 
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FIG. 11. The melting line T^(B), for e = 0.02,0.1,0.5, from 
the Lindeman criterion, assuming a temperature independent 
penetration depth A. At high fields, the melting tempera- 
ture decreases with the magnetic induction as B""'^"^, as in- 
dicated by the asymptotic dashed lines, and linearly with the 
anisotropy parameter e. The dotted-dashed lines show the 
melting lines resulting from the local elastic moduli. While 
the latter agrees quite well with the reentrant low field part of 
the curve, it yields a very different behavior at higher fields. 



In Figs. ^ |T^, and |l^, we depict the resulting melt- 
ing lines for different values of the anisotropy parameter 
e = 0.02 (BiSCCO), and e = 0.1 (YBCO). Note that 
some of the data are presented in a log-log or in a log- 
linear plot in order to emphasize the scaling with the 
magnetic field in the high field region, as well as to bet- 
ter visualize the lower part of the melting line, i.e. the 
reentrant portion of the phase diagram. The dotted line 
in the last two figures represents the upper critical field 
' ''{T) = (t)o{l - T/Tc)/2n^^{0), with = 93°K and 



ttMF 



f(0) = UA. Close to H^'^iT) the amplitude of the or- 
der parameter is strongly reduced due to the overlap of 
the vortex cores and these results are no longer valid. 
However, well below this line, the London limit provides 
a good description of the system. In Fig. nil we repre- 
sent the results obtained assuming A(r) — A(Oj — 1400A, 
for e = 0.02,0.1,0.5. At high fields a/A < 1, the melt- 
ing temperature decreases as jjj a,ll cases, consis- 
tent with the prediction oLJL—!^ _B^^/^, common to the 
above mentioned modela23Gil'E3. We have also checked 
that the melting temperature decreases linearly with e, 
for various values of a/ A within the high field region, 
i.e. Tm ^ eB~^-^^. This is another feature which is in 
agreement with previous predictions, and with available 
experimental data. In this figure we also encounter some 
unreasonably high melting temperatures which are the 
consequence of the assumption A(T) = A(0). 
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FIG. 12. The melting line r,„(B), for e = 0.02,0.1, assuming 
a mean-field temperature dependent penetration depth A. The 
dashed lines show the results obtained assuming A = A(0) up 
to the point in which they are still reasonable. 

Figure |l^ shows the results obtained for the melting 
line r„ assuming A(r) = A(0)(1 - T/Tc)-^/^. We have 
used the values of Tc = 90 for BiSCCO (e = 0.02), and 
Tc = 93 for YBCO (e = 0.1). As in the previous fig- 
ure, within the high field region of the phase diagram, 
the melting temperature decreases as the external field 
is increased. Quantitative and qualitative differences be- 
tween the new curves and the corresponding results for 
A(r) = A(0) (depicted as dashed fines in the figure) can 
be observed even at low temperatures. Near the critical 
temperature, the low-field portion of the melting curves 
has a negative slope (as for their high field counterparts), 
but at low temperatures the slope is positive, so that the 
melting temperature decreases with the external field at 
very low densities of flux lines. 

At the scale of Fig. |l^, the high and low field portions 
of the melting curve Tm{B) appear almost discontinuous 
close to Tc- However, as indicated in the blow-up of this 
region in Fig. |l^, this is in reality a sharpj-but, contin- 
uous turn around. Following the literaturecJ~Ea, we fit 
the melting fields i?m(T) to power laws in the reduced 
temperature t = 1 — TjTc, with an exponent of a. The 



dashed lines in Figure 13 indicate the extent of this power 
law regime for each value of e. We note that even though 
we are using the input value of Tc, the power law regime 
extends at most over two decades. We naturally obtain 
different exponents for the upper and lower branches of 
the melting curve. The apparent exponent of the upper 
branch actually depends on the parameter e, taking the 
value of a+{e = 0.1) = 1.91 for 'YBCO', and a slightly 
smaller value of a+{e = 0.02) = 1.67 for 'BiSCCO'. Both 
curves overlap in the reentrant region at very dilute con- 
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centrations where we can fit to an effective exponent of 
a- = 0.75. 

The phase diagram of Fig. |l^ is quahtatively similar to 
the predictions of Ref. ^, where the power law forms for 
the phase boundary close to Tc were also first discussed. 
In particular, it is reasonably straightforward to estimate 
the shape of the phase boundary in the low field region, 
where the leading contributions to the transversal modes 
At come from the shear modulus cge (which decays ex- 
ponentially with a/ A, and does not depend on e), and 
from the last term in Eq. (^) . As expected for very low 
densities of flux lines, the last two terms in this equa- 
tion tend to the single line tension £{kz) (see Eq. (|T9|)) 
which, for small values of fc^, is in turn dominated by 
the e-independent magnetic contribution. This leads to 
i3~ ~ A(T)~^ ln~^(A(T)), and consequentlypto an expo- 
nent of a_ = 1 with logarithmic correctionsB. The effec- 
tive exponent of a_ « 0.75 obtained in Fig. ^is different 
from the a_ = 1 expected on the basis of our mean-field 
form for A(r) indicating the difficulty of determining the 
true exponent from such fits. 
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FIG. 13. The low-field region of the reentrant melting line 
plotted as a function of (1 - T/TJ, for e = 0.02,0.1. 
Close to the critical temperature Tc, the melting induction 
grows as a power of the reduced temperature 1 — T/Te, as 
indicated by the asymptotic dashed lines. As in Fig. hi], the 
dotted-dashed lines show the melting lines resulting from the 
local elastic moduli. 

For the upper portion of the mean-field line, the melt- 
ing field is estimatedcl to scale as ~ A(r)~^, leading 
to a+ = 2 for our choice of a mean- field A(T). While the 
observed value of Q!+(e = 0.1) = 1.91 is not far from this 
expectation, the smaller value of a+(e = 0.02) = 1.67 in- 
dicates the difficulty of finding a good asymptotic regime. 
Naturally, in close vicinity of the critical point we can 



no longer use the mean field forms for the divergence of 
A (X ^ oc Using the scaling forms, ^ (x \t\~^ oc A^, 

with the XY critical exponent of v w 2/3, leads to 
a+ = 2a_ = 2i>El. Indeed, the values of the exponent 
of the upper branch (a+ ^ 1.33 — 1.36) reported from 
measurenaents of the melting transition in single crystals 
of YBCOE£l"Ea, are consistent with this prediction. How- 
ever, given the difficulties of extracting exactly known 
exponents from the data of Fig. |l3|, we may well ques- 
tion the robustness of this procedure. 

We would like to emphasize that most of these results 
cannot be obtained fr om t he local values of the elastic 
moduli given in Sec. IVB. The results obtained from 



such a local elastic approximation are plotted as dot- 
dashed lines in Figs. |l^ and ^ The latter curves overlap 
in the high field region, where neither the local tilt, nor 
the shear modulus depend on e. Rather than decreas- 
ing as r„i ~ B~^l'^ ^ the resulting melting temperature 
reaches a constant value for high _B, i.e. Tm ~ -B"e". On 
the other hand, the local melting temperature provides 
a very good approximation to the melting line at very 
dilute concentrations, in the reentrant low field portion 
of the curve. 



VII. CONCLUSIONS 



In conclusion, we have obtained the elastic moduli of 
the flux line lattice in a systematic and detailed way 
which is not restricted to the most commonly consid- 
ered continuum limit. The transversal and longitudi- 
nal harmonic eigenvalues have been computed numeri- 
cally for different areal densities of flux lines, and as a 
function of Q (within the irreducible Brillouin zone) and 
kz- Several features emerge from the analysis of our re- 
sults: (i) There is a weak angular dependence of the 
harmonic eigenvalues which becomes more pronounced 
on approaching the BZ boundary, (ii) Throughout the 
BZ, transversal modes are less costly than longitudinal 
modes, and are the main cause of lattice fluctuations, 
(iii) Not surprisingly, both eigenvalues increase with \kz , 
and modes with higher fcz are more costly, (iv) Rather 
surprisingly, due to a rapid decrease of C44, the energy of 
a transversal mode with a nonzero fc^ actually goes down 
with Q. The minimum cost occurs for a finite Q which 
depends on Afc^, and the density of flux lines, (iv) For a 
large portion of the IBZ, both eigenvalues exhibit a linear 
dependence on (aQ)^, whose slope depends on the local 
shear modulus cqq which can be calculated as a function 
of a/A. 

Some of the above features are qualitatively well ac- 
counted for by the nonlocal continuum limit results recal- 
culated in Sec. IV. Nevertheless, beyond a characteristic 
value of Afcz, there are major differences between this an- 
alytic form and those calculated numerically. Guided by 
our results, we propose analytic corrections to the non- 
local continuum limit, which fit quite well the behavior 
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of the elastic eigenvalues throughout the London regime. 
Other limiting forms for the elastic energy often quoted 
in the literature, such as the single line, local, and lo- 
cal continuum approximations, are reconsidered within 
our framework, and their range of validity is examined. 
We have plotted some of these limits, together with our 
numerical results, for better ease of comparison. 

We also consider different values of the anisotropy pa- 
rameter e, in particular corresponding to values quoted 
in the literature for high temperature superconductors 
such as YBCO or BiSCCO. The latter is the paradigm 
of a highly anisotropic superconductor with softer elastic 
eigenvalues. This fact can be corroborated in our analysis 
and is ultimately responsible for strengthening the mag- 
nitude of thermal fluctuations. We have made use of our 
proposed analytical expression to calculate the extent of 
thermal fluctuations for different e. The full form of the 
melting curve is then obtained using the Lindeman cri- 
terion, capturing the following salient features of T,„(_B): 
(i) It decreases at high temperatures with the magnetic 
field, approximately as B~'^-^^. (ii) It decreases with the 
anisotropy parameter, and consequently the molten vor- 
tex liquid covers a large fraction of the equilibrium phase 
diagram at small e. (iii) There is reentrant melting at low 
fields due to the weak interactions between widely sepa- 
rated flux lines. The reentrant phase boundary is itself 
a non-monotonic function of temperature, (iv) Close to 
the critical temperature Tc, both branches of the melting 
line Bm can be fitted to power laws in the reduced tem- 
perature 1 — T/Tc- However, the power law is observed 
at most over two decades, and the resulting effective ex- 
ponents arc different from the expectations from mean 
field theor}£_casting doubt upon the effectiveness of this 
methodfeaia. 

Another interesting potential extension of this work is 
to include entropic contributions to the free energy, in or- 
der to explore the fluctuation induced effects which have 
been proposed to be .important in the low field region 
of the phase diagramEZl. Our analysis has been applied 
to a conventional s-wave superconductor with a trian- 
gular lattice of flux lines. A square lattice of flux lines 
is also possible in d-wave superconductors, where simi- 
lar approach and phenomenology can be carried out, al- 
though the details of the harmonic energy, as dictated 
by symmetry, will be different. Again, calculation of the 
entropic contributions to the free energy may provide a 
better estimatp-pf the transition between triangular and 
square latticesH. 
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